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ABSTRACT 

The redshifted 21-cm signal from neutral hydrogen (HI) is an important fu- 



j> . ture probe of the the high redshift universe. We have analyzed 610 MHz 

■ GMRT observations towards detecting this signal from z — 1.32. The multi- 

frequency angular power spectrum Ci(Av) is used to characterize the statis- 
tical properties of the background radiation across angular scales ~ 20 to 



10 , and a frequency bandwidth of 7.5 MHz with resolution 125 kHz. The 
measured Ce(Av) which ranges from 7 mK 2 to 18 mK 2 is dominated by fore- 
grounds, the expected HI signal C™(Ai/) ~ 10~ 6 — 10~ 7 mK 2 is several 
orders of magnitude smaller and detecting this is a big challenge. The fore- 
grounds, believed to originate from continuum sources, is expected to vary 
smoothly with Av whereas the HI signal dccorrclates within ~ 0.5 MHz and 
this holds the promise of separating the two. For each £, we use the inter- 
val 0.5 < Av < 7.5 MHz to fit a fourth order polynomial which is subtracted 
from the measured Ce(Au) to remove any smoothly varying component across 
the entire bandwidth Av < 7.5 MHz. The residual Ce(Ais), we find, has an 
oscillatory pattern with amplitude and period respectively ~ 0.1 mK 2 and 
Av = 3 MHz at the smallest I value of 1476, and the amplitude and period 
decreasing with increasing I. Applying a suitably chosen high pass filter, we are 
able to remove the residual oscillatory pattern for I — 1476 where the residual 
Ci(Av) is now consistent with zero at the 3er noise level. Based on this we 
conclude that we have successfully removed the foregrounds at t = 1476 and 
the residuals are consistent with noise. We use this to place an upper limit on 
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the HI signal whose amplitude is determined by xmb (C™(Aj/) oc [xhi&] 2 ), 
where xhi and b are the HI neutral fraction and the HI bias respectively. A 
value of xmb greater than 7.95 would have been detected in our observation, 
and is therefore ruled out at the 3cr level. For comparison, studies of quasar 
absorption spectra indicate xhi ~ 2.5 x 10~ 2 which is ~ 330 times smaller than 
our upper limit. We have not succeeded in completely removing the residual 
oscillatory pattern, whose cause is presently unknown to us, for the larger i 
values. 

Key words: cosmology: observations, cosmology: diffuse radiation, cosmol- 
ogy: large-scale structures 

1 INTRODUCTION 

Detecting redshifted 21-cm radiation from neutral hydrogen (HI) at high redshifts is of 
considerable interest in cosmology. At redshifts z < 6, the bulk of the neutral gas is in 
clouds that have HI column densities in excess of 2 x 10 20 atoms/cm 2 (Lanzetta et al. 1995; 
Storrie-Lombardi, McMahon & Irwin 1996; Peroux et al. 2003). These high column density 
clouds are observed as damped Lyman-o; absorption lines seen in quasar spectra. The analysis 
of quasar spectra indicate that the ratio of the density p ga , s (z) of neutral gas to the present 
critical density p co of the universe has a nearly constant value Q gas (z) ~ p gas (z)/ 'p co ~ 10 -3 , 
over a large redshift range 0.5 < z < 5.0 (Storrie-Lombardi, McMahon & Irwin 1996; 
Rao & Turnshek 2000; Peroux et al. 2003; Prochaska & Herbert-Fort 2004; Rao, Turnshek 
& Nestor 2006; Kanekar et al. 2009). The redshifted 21cm radiation from the HI in this 
redshift range will be seen in emission. The emission from individual clouds (< 10 pJy) is too 
weak to be detected with existing instruments unless the image is significantly magnified 
by gravitational lensing (Saini, Bharadwaj & Sethi 2001). The collective emission from 
the undetected clouds is present as a very faint background in all radio observations at 
frequencies below 1420 MHz. The fluctuations in this background radiation carry an imprint 
of the HI distribution at the redshift z where the radiation originated. The possibility of 
detecting this holds the potential of providing us with a new observational probe of large- 
scale structures (Kumar, Padmanabhan & Subramanian 1995; Bagla, Nath & Padmanabhan 
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1997; Bharadwaj, Nath & Sethi 2001; Bharadwaj & Sethi 2001; Bharadwaj & Pandey 2003; 
Bharadwaj & Srikant 2004; Ali,Bharadwaj & Pandey 2005; Wyithe & Loeb 2008, 2007; 
Bagla & Khandai 2009; Khandai, Datta & Bagla 2009). In a recent paper Pen et al. (2009) 
report a detection of the post-reionization HI signal through the cross-correlation between 
the HIPASS and the 6dfGRS data. 

Observations of redshifted 21-cm radiation can in principle be carried out over a large 
redshift range starting from the cosmological Dark Ages through the Epoch of Reionization 
to the present epoch (Bharadwaj and Ali 2005), allowing us to trace out both the evolution 
history of neutral hydrogen as well as the growth of structures in the universe. Redshifted 
21-cm observations also hold the potential of allowing us to probe the expansion history 
of the universe (Mcquinn et al. 2006; Chang et al. 2008; Visbal, Loeb & Wyithe 2008; 
Bharadwaj, Sethi & Saini 2009). 

The Giant Meter Wave Radio Telescope (GMRT 1 ; Swarup et al. 1991), currently oper- 
ating at several frequency bands in the frequency range 150 to 1420 MHz is well suited for 
carrying out observations towards detecting the HI signal over a large redshift range from 
z ~ to z ~ 8.3 and angular scales of ~ 10" to ~ 1°. In this paper we report results from 
the analysis of 610 MHz observations towards detecting the redshifted 21-cm signal from the 
cosmological HI distribution at z — 1.32. 

We have characterized, possibly for the first time, the statistical properties of the back- 
ground radiation at 610 MHz across ~ 20" to 10' angular scales and a frequency bandwidth 
of 7.5 MHz with a resolution of 125 kHz using the multi-frequency angular power spectrum 
Ce(Av) (hereafter MAPS; Datta, Roy Choudhury & Bharadwaj 2007). This jointly charac- 
terizes the angular (£) and frequency (Az/) dependence of the fluctuations in the 610 MHz 
radiation in the field of view of our observation. Foregrounds from different astrophysical 
sources are expected to be a few orders of magnitude larger than the predicted 21-cm signal 
(Shaver et al. 1999; Di Matteo et al. 2002; Oh and Mack 2003; Santos et al. 2005; Wang 
et al. 2006; Ali, Bharadwaj & Chengalur 2008) and our 610 MHz GMRT observations are 
expected to be nearly entirely dominated by foregrounds which are predicted to be at least 
a thousand times larger than the HI signal. Separating the HI signal from foregrounds is the 
most important challenge for cosmological redshifted 21-cm observations. 

The foregrounds are believed to have a smooth continuum spectra and the contribution to 
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Ce(Av) is expected to vary very slowly with Av across the band (7.5 MHz) of our analysis. 
The contribution from the HI signal decorrelates very rapidly with increasing Av and is 
expected to be uncorrected beyond Av = 0.5 MHz at the angular scales (£ = 10 3 to £ = 
3 x 10 4 ) of our analysis. This property of the signal holds the promise of allowing us to 
separate the signal from the foregrounds. In this paper we propose and implement a technique 
that uses polynomial fitting in Av to subtract out any smoothly varying component from 
the measured Cp(Av). The residuals are expected to contain only the HI signal and noise. 
The target of the present work is to test if the polynomial subtraction successfully removes 
the foregrounds to a level such that the residuals are consistent with noise. The noise in the 
current observation is considerably larger than the Hl-signal and longer observations would 
be needed for detecting the HI signal. 

The present work closely follows an earlier paper (Ali, Bharadwaj & Chengalur 2008) 
which analyzed 150 MHz GMRT observations. We note that the prospect of detecting the 
redshifted 21-cm signal considerably increases at higher frequencies (eg. 610 MHz) where the 
foreground contribution and noise are both smaller. Further, the problem of man made radio 
frequency interference is considerably more severe at 150 MHz as compared to 610 MHz. 

A brief outline of the paper follows. Section 2 describes the observation and data analysis, 
Section 3 presents the visibility correlation technique that we use to estimate Ce(Av) and 
also presents the estimated values, Sections 4 and 5 present model predictions for the HI 
signal and foregrounds respectively, while Section 6 describes our proposed technique of 
foreground removal and finally Section 7 contains results and conclusions. 

2 GMRT OBSERVATIONS AND DATA ANALYSIS 

The GMRT has 30 fixed antennas each of diameter 45 m. 14 of which are randomly dis- 
tributed in a central square 1.1km x 1.1km in extent, while the rest of the antennas are 
distributed approximately in a 'Y' shaped configuration. The shortest antenna separation 
(baseline) is around 60 m including projection effects while the largest separation can be as 
long as 26 km. The hybrid configuration of the GMRT gives reasonably good sensitivity to 
probe both compact and extended sources. 

The observed field of view is centered on 0:2000 — 12 ft 36 m 49 s , #2000 = 62°17'57" which is 
situated near Hubble Deep Field North (HDF-N) (a 2 ooo = 12 fc 36 m 49.4 a , <5 20 oo = 62°12'58" ). 
The galactic co-ordinates of the observed field is I = 125.87°, b = 54.74°. The sky temper- 
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ature determined at this location is 20 K in the 408 MHz Haslam et al. (1982) map. The 
observations were carried out over three days from 4 th to 7 th September, 2002 and the total 
observation time was almost 30 hours (including calibration). The observation had a center 
frequency of 618 MHz, and a total band-width of 16 MHz divided into 128 frequency chan- 
nels, each 125 kHz wide. The integration time was 16 seconds and visibilities were recorded 
for two orthogonal circular polarizations. The calibrator source 3C147 and 3C286 were used 
for flux calibration and 1313+675 was used for phase calibration. The phase calibrator was 
observed every half hour to correct for temporal variations in the system gain. We have used 
the Astronomical Image Processing Software (AIPS) to analyze the recorded visibility data. 
The flux of these two flux calibrator was estimated by extending the Baars scale (Baars et. 
al. 2005) to low frequencies using the AIPS task 'SETJY'. Standard AIPS tasks were used 
to flag all data that could be visually identified as being bad. The entire lower sideband data 
was found to be bad and was discarded from the subsequent analysis. Data from different 
days were calibrated and flagged separately and then combined using the AIPS task 'DB- 
CON'. We find that the channels near the edge of the band are relatively noisy and hence 
only the 100 central channels were used in the subsequent analysis. 

An initial 2D image of the field of view (FOV) showed four bright sources with consider- 
able imaging artifacts. To improve our image quality , initially we have subtracted out the 
clean components (CC) of these bright sources by moving them to the phase center using 
appropriate RA-SHIFT and DEC-SHIFT within AIPS. Then, we add back the brightest 
source and use this for three rounds of self-calibration with time intervals of 3 and 2 minutes 
for phase calibration and finally 20 minutes for amplitude and phase calibration and then 
subtract out the brightest source again. The same process is followed for rest of the bright 
sources. Subsequent to this, we have also subtracted out all the weaker sources from our 
FOV and used the AIPS task 'TVFLG' to flag out any bad visibility. Finally, we have col- 
lapsed all frequency channels and clipped the resulting visibilities at 0.07 Jy. At each stage 
the same calibration and flag tables were also applied to the original 100 channel data which 
contains all the sources. 

The large field of view (# FW hm = 43') of the GMRT at 610 MHz leads to considerable 
errors if the non-planar nature of the GMRT antenna distribution is not taken into ac- 
count. We use the three dimensional (3D) imaging feature ( Perley 1999) in the AIPS task 
'IMAGR' in which the entire field of view is divided into multiple sub-fields, each of which 
is imaged separately. Here, a 1.5° x 1.5° FOV was imaged using 163 facets. The presence 
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Figure 1. Our continuum image of bandwidth 12.5 MHz centered at 618 MHz. The FOV was imaged using 163 facets which 
have been combined using the AIPS task FLATN. The rms off-source noise is 60/^Jy/Beam. Most of the extended features 
visible in the image are imaging artifacts around the four bright sources. 



of a large number of sources in the field allows us to do self calibration loops to improve 
the image quality. We have applied 4 rounds of self calibration, the first three only for the 
phase and the final round for both amplitude and phase. The time interval for the gain 
correction was chosen as 3, 2, 1 and 20 minutes for the successive self calibration loops. At 
every stage the calibration tables were applied to the original 100 channel data. The 100 
channels were collapsed into 10 channels which were used to make a continuum image of 
the entire FOV. Figure 1 shows our continuum image of band width 12.5 MHz centered at 
617 MHz. The synthesized beam has FWHM 8.2 x 5.8 . The off source rms. noise in the 
image was ~ 60 /xJy/Beam and the image quality had improved considerably. The rms noise 
around the brighter sources is higher, using the AIPS task 'TVSTAT' we notice that it is 
around 100 /iJy/Beam. The brighter sources are also found to be accompanied by a region 
of negative flux density, these are presumably the results of residual phase errors which were 
not corrected in our self calibration process. The maximum and minimum flux density in 
the final image are ~ 250mJy/Beam and ~ — 2.8mJy/Beam respectively. 

The subsequent analysis was done using the calibrated visibilities of the original 100 
channel data with all the sources. The data contains 510528 baselines, each of which has 
visibilities for 2 circular polarizations. The baselines are in the range 200 A to 20 kA. The 
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Figure 2. Histogram showing the distribution of the real parts of all the 33233698 measured visibilities. The same is shown 
on a linear scale (left panel) and a log-linear scale (right panel). A Gaussian with mean 0.683 mjy and rms 0.251Jy, values 
calculated for the real parts of the measured visibilities, is plotted as a dotted line in both panels. While the Gaussian fits the 
data very well at small amplitudes, there is a discrepancy at high amplitudes (~ 0.75Jy) which is visible only in the right panel. 

visibilities from the two polarizations were combined (V = [Vrr + Vll]/2) for the subse- 
quent analysis. The final calibrated data contains 33233698 visibilities. The real part of the 
visibilities have a mean of 0.68 mjy and rms of 0.25 Jy. Similarly, the imaginary part has a 
mean of 2.34 mjy and rms of 0.25 Jy. 

It is often convenient to assume that the visibilities have a Gaussian distribution. The 
distribution of the real part of the visibilities is shown in a histogram (Figure 2). We find 
that a Gaussian gives a reasonably good fit to the data within 3cr, which contains the bulk 
of the data. The number counts predicted by the Gaussian falls much faster than the data 
at large visibility values | Re(V) |> 0.75 Jy. The imaginary part of the visibilities show a 
similar behavior. 

Deviation from Gaussian statistics is expected to mainly effect the error estimate in the 
visibility correlation. We expect this effect to be small , since the discrepancy is for only 
small fraction of visibilities. 



3 VISIBILITY CORRELATIONS AND THE ANGULAR POWER 
SPECTRUM 

The aim here is to quantify the statistical properties, in angle and frequency, of the 610 MHz 
sky signal. For a frequency z/, the angular dependence of the brightness temperature distri- 
bution on the sky T(u, n) may be expanded in spherical harmonics as 

T(z/,n) = ]T a tm (u) Y tm (n) . (1) 

£,m 
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The multi-frequency angular power spectrum (MAPS; Datta, Roy Choudhury & Bharadwaj 
2007), which jointly characterizes the dependence on angular scale and frequency separation, 
is defined as 

Q(Av) = Q(v, v + Av) = (a lm (v) a\ m {y + Av)) . (2) 

Here i refers to the angular modes on the sky. The sky signal is assumed to be statistically 
isotropic. We also assume that for the relatively small bandwidth of our observation (Av <C 
v), the frequency dependence can be entirely characterized through Av whereby we do not 
explicitly show v as an argument in eq. (2). 

We use the correlation between pairs of visibilities V(U, v) and V(U + AU, v + Av) 

V 2 (U, v; U + AU, v + Av) ee (V(U, v)V*(U + AU, v + Av)) (3) 

to estimate C e (Av). The correlation of a visibility with itself is excluded to avoid a positive 
noise bias in the estimator. Ali, Bharadwaj & Chengalur (2008) as well as Dutta et al. (2009) 
contain detailed discussions of the estimator and we highlight only a few salient features here. 

The GMRT primary beam pattern is well approximated by a Gaussian A{6) = e~ e2 / e o 
where 8 « 0.6 x #fwhm = 25'. 8 (6> F whm = 43') at 610 MHz. In a situation where AU 
is small such that | AU |< (ttOq)" 1 = 42.4 A (9 in radians), the expected correlation 
V 2 (U, v; U + AU, v + Av) in eq. (3) does not depend on AU whereby we may express it as 
V 2 {U,Av). Further, if U > (vr^)- 1 we have 

V 2 (U,Av) = ^ (^j C 2 ^{Av)Q{Av). (4) 

where B v = 2v 2 k B T/c 2 is the specific intensity of black-body radiation in the Rayleigh- Jeans 
limit. Both #0 and(||r) depend on the frequency. In our analysis we treat these as constants 
with the value being evaluated at 610 MHz. It is possible to incorporate the effect of the 
Av dependence of 6 and^^f ) through the function Q(Av) in eq. (4). This is expected to 
be a slowly varying function of Av with a variation of ~ 1% across the Av range of our 
observation. We have not explicitly considered the function Q(Av) in our present analysis. 
This is expected to introduce an extra, slowly varying Av dependence in the estimated 
Ci(Av). This slowly varying Av dependence, as we shall discuss later, can be included 
in the foreground model and separated from the HI signal which varies rapidly with Av. 
Equation (5) gives the final expression that we use to estimate the angular power spectrum 
(MAPS) 
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Figure 3. The real (upper) and imaginary (lower) parts of the observed visibility correlation V2(U, 0). The real part of V2(U, 0) 
may be interpreted as Ci where £ = 2n U (labeled on the right and top) . The 1 — a error-bars shown have contributions from 
both the cosmic variance and system noise. 



C 2nU (Au) = 87 {— x V 2 (U, Av) . (5) 

In our analysis we have correlated only pairs of baselines which satisfy the condition 
| AU |< 10A. We have restricted the analysis to baselines 200 A < U < 5000 A. To test if 
the visibility correlation is actually independent of AU we have also considered | AU |< 5A 
and 20 A. The results are unchanged for 20 A and they are rather noisy for 5A, there being 
very few baseline pairs to correlate. 

The measured V 2 (U, Au) will, in general, have real and imaginary parts (Figure 3). As 
seen in eq. (5), the expectation value is predicted to be real, the expectation value of the 
imaginary part being zero. We use the real part of the measured V 2 (U, Au) to estimate 
Ce(Au) through eq. (5). A small imaginary part arises due to the noise in the individual 
visibilities. This introduces random fluctuations in both the real and imaginary parts of 
the measured V 2 (U,Au). Figure 3 shows the measured V 2 (Au) and the inferred Ce(Au) for 
Av = 0. As expected, the imaginary part is much smaller than the real part of V 2 (U, Au). 
Note that we use the notation C e = Ce(Au = 0). 

We next consider the expected statistical fluctuations (error) in V 2 (U, Av). The total error 
has two parts i.e., system noise and the cosmic variance. The total variance [AV 2 (Ui, Az/)] 2 
can be calculated as 



10 A. Ghosh, S. Bharadwaj, S. S. Alt and J. N. Chengalur 



,, (ivy , [v 2 (u„M? 



[AViW.A^^ + L^-^ (6) 

where (A^ 2 ) = (NN*) is the variance of the noise contribution N in the visibilities that 
we use in our analysis, N P is the total number of baseline pairs that contribute to V 2 (U, Au) 
and N E is the number of independent estimates of V^U, Au). Here (A^ 2 ) = a 2 where a is 
the rms. noise, for a single polarization, in the real part (or equivalently the imaginary part) 
of a visibility. The value of a is expected to be (Thompson, Moran & Swenson 1986) 

\/2kBT S y S . . 

a = / ( ' ) 

A eff VAuAt 

where T sys is the total system temperature, ks is the Boltzmann constant, A e ff is the 
effective collecting area of each antenna, Av is the channel width and At is correlator 
integration time. For the parameters of our observations, T sys 100 K, TT sys k B i 'A e ff = 
300 Jy, Au = 0.125 MHz and At = 16 s we have a 2 = 2.25 x 10~ 2 Jy 2 . In our analysis we 
have used (iV 2 ) = 1.25 x 10 _1 Jy 2 which is the sum of the variance of the real and imaginary 
components of the measured visibilities. In our observation the total error is dominated by 
the cosmic variance which is a few orders of magnitude larger than the system noise in the 
entire U range that we have considered. 

The Au dependence of C e (Au) is shown in Figure 4. We have considered U values below 
1.25 kA where £ = 2nU . As discussed later, the HI signal falls at U > 1 kA which is why 
we have not considered baselines larger than 1.25 kA. We find that for nearly all the values 
of £ shown in the figure the variation in Cg(Av) with Au is roughly between 0.2 mK 2 to 
0.6 mK 2 across the 7.5 MHz band. The fractional variation in C e (Au) ranges from 1.5 % to 
3.6 %. We note that an oscillatory pattern is visible in Ce(Au) at nearly all values of I. The 
pattern is most pronounced at the lower £ values. The error-bars shown in Figure 4 include 
only the system noise contribution. The measured Cg(Au) is expected to be dominated by 
foregrounds which are believed to be largely independent of Au. For a fixed £ the cosmic 
variance then is expected to introduce the same error (independent of Au) across the entire 
band. As a consequence we do not consider the cosmic variance for the Au dependence 
shown in Figure 4. 

The two dimensional (2D) Fourier transform relation between the sky brightness and 
the visibilities assumed in deriving eq. (4) is not strictly valid for GMRT's field of view 
(#fwhm — 43 ). In addition to uv which are the components of the baseline in the plane 
normal to the direction of observation, it is also necessary to consider w the component along 
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Figure 4. The measured C't(Av) as a function of Au for the £ values shown in each panel. The error-bars include only 3<r 
system noise. The solid line shows the 4 th order polynomial fits. The details of the polynomial fitting procedure are discussed 
in Section 6. 



the observing direction. To assess the impact of the w term we have repeated the analysis 
using only a limited range of baselines for which w < 0.5 x U. We find that limiting the 
maximum w value does not make any qualitative change in our results. 



4 THE EXPECTED REDSHIFTED HI 21 CM SIGNAL 

Our observing frequency v — 610 MHz corresponds to a redshift of z =1.32 for the HI 
21-cm radiation. Observations of Lyman-a absorption lines seen in quasar spectra indicate 
that the ratio of the density p gas (z) of neutral gas to the present critical density p crit of 
the universe has a nearly constant value p ga , s (z) / p crit ~ 10 -3 , over a large redshift range 
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< z < 3.5. This implies that the mean neutral fraction of the hydrogen gas is xhi = 
50 Q gas h 2 (0.02/Qi ) h 2 ) = 2.45 x 10~ 2 which we adopt for our analysis. The redshifted 21cm 
radiation from the HI will be seen in emission as a very faint background in our observation. 
The fluctuations in this background with angle and frequency is a direct probe of the HI 
distribution at the redshift z— 1.32 where the radiation originated. We calculate the MAPS 
for the redshifted 21-cm signal (Datta, Roy Choudhury & Bharadwaj 2007) using 

d(Au) = — / dfc|| cos(fc|| r' v Au) P m (k) , (8) 
nrf, Jo 

where the three dimensional wave vector k has been decomposed into components k\\ and 
l/r u , along the line of sight and in the plane of the sky respectively. The comoving distance 
r v is the distance at which the HI radiation originated. Note that (1 + z)^ 1 r v = dpj^z) is the 
angular diameter distance and r' v = dr v /dv. The temperature occurring in eq. (8) is given 
by 

f(z)=40mK (! + ,)• fB^VH^, (9) 



0.02 J V h J H(z) ' 

and -Pffl(k) is the three dimensional power spectrum of the "21 cm radiation efficiency in 
redshift space" (Bharadwaj and Ali 2005) which in this situation is given by 

P m (k)=x 2 m b 2 (l + Pfi 2 ) 2 p(k). (10) 

The term (1 + (3/j 2 ) 2 arises because of the HI peculiar velocities (Bharadwaj, Nath & Sethi 
2001; Bharadwaj & Ali 2004), which we assume to be determined by the dark matter. This 
is the familiar redshift space distortion seen in galaxy redshift surveys, where j3 is the linear 
distortion parameter and \x = k\\/k. On the large scales of interest here, it is reasonable 
to assume that HI traces the dark matter with a possible linear bias b, whereby the three 
dimensional HI power spectrum is b 2 P(k), where P(k) is the dark matter power spectrum at 
the redshift where the HI signal originated. Unless mentioned otherwise, we use the values 
(fi m o, ^ao, h, <7 8 , n s ) = (0.3,0.7,0.7,1.0,1.0) for the cosmological parameters and b = 1 for 
the bias throughout. 

Figure 5 shows Ce = Ci(Av = 0) for the expected HI signal. The Av dependence has 
been shown (Figure 6) through the frequency decorrelation function ki(Av) which is defined 

as 

, . x CAAv) , . 

This function quantifies how quickly the HI signal decorrelates as we increase the frequency 
separation Av, with the signal being correlated and uncorrelated when ne(Av) ~ 1 and 
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Figure 5. The predicted Ct for the 610 MHz redshifted 21-cm signal. 



Ki(Au) ~ respectively. Bharadwaj & Pandey (2003) have proposed an analytic formula 
that approximates the visibility correlation V 2 (U,Av) for the GMRT 610 MHz HI signal. 
Using this we obtain the analytic expression 

_ /A , , /1000y / AzA . /2tt AzA 

with A = 8.0 x 10- 7 mK 2 , 7 = 1.2, 6 £ = 0.48 x (1000/f) a8 MHz and d e = 1.8 x 
(1000/f) L2 MHz which matches the numerically calculated values of Ce(Av) at around 
~ 10% for the £ and Av range of our interest. 

We find that the predicted Ci declines rapidly (oc i^ 1 ' 2 ) with increasing t. Based on this 
we have restricted our analysis to the I range 1000 to 5000 where the signal is expected to 
be strongest. The signal would be larger at t < 1000, but the GMRT's field of view restricts 
us from accessing these £ values. The measured C e values are around 10 7 to 10 s times larger 
than the predicted HI signal. The predicted signal decorrelates rapidly with increasing Av 
and Kt(Av) falls by 90% or more (Ke(Au) < 0.1) at Av = 0.5 MHz. The value of Av where 
Ke(Av) falls by 90% is smaller for larger values of i. Further, we find that the expected 
HI signal is anti-correlated at large values of Av (~ 0.8 MHz) where Ke(Au) has a small 
negative value. 
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Figure 6. The predicted K((Av) for the 610 MHz 21-cm signal at three different £ values in the range 1000 to 5000. 



5 FOREGROUND MODEL PREDICTIONS 

The radiation coming from different astrophysical sources, other than the HI signal, con- 
tribute to the foreground radiation. Here we mainly focus on the two most dominant fore- 
ground components namely extragalactic point sources and the diffuse synchrotron radiation 
from our own Galaxy. The free-free emissions from our Galaxy and external galaxies (Shaver 
et al. 1999) make much smaller contributions though each of these is individually larger than 
the HI signal. We have modeled the MAPS for each foreground component as 



where A, 7 and k^(Az/) are the amplitude, the power law index and frequency decorrelation 
function respectively. The different foreground components considered here are all continuum 
radiation which are known to vary smoothly with frequency. For each component we denote 
the spectral index using a, whereby the amplitude scales as A oc u 2a . The value of A, 
whenever used in this paper, is at a fixed frequency of 610 MHz. The continuum nature of 
the foreground components also implies that we expect fi^(Az/) to be of order unity and vary 
smoothly with Av. The foregrounds will remain correlated across the frequency band of our 
observation, unlike the HI signal which decorrelates rapidly within Au = 0.5 MHz. Given 
the absence of any direct observational constraints on k^{Av) for any of the foreground 
components at the angular scales and frequencies of our interest, we do not attempt to 




(13) 
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make any model predictions for this quantity beyond assuming that it varies smoothly with 
Au across the frequency band of observation. In the subsequent discussion we focus on 
model predictions of A and 7 which are tabulated in Table 1 for the different foreground 
components. 

Extra-galactic point sources are expected to dominate the sky at 610 MHz. We have 
estimated the point source contribution using the 610 MHz differential source count from 
Garn, Green, Riley et al. (2008). This is the average differential source count of 610 MHz 
GMRT observations in three different fields of view namely the Spitzer extragalactic First 
Look, ELAIS-N1 and Lockman Hole surveys. The differential source count in the flux range 
of their observation (~ 0.3 to 200 mJy) is well fitted by a single power law 
(W 1259 / S V 1 ' 84 

~dS ~ JySr ' \Uy) ' ( ' 

We have assumed that the same power law also holds for the fainter sources below the 

detection limit. 

Point sources make two distinct contributions to the angular power spectrum, the first 
being the Poisson noise due to the discrete nature of the sources and the second arising 
from the clustering of the sources. The Poisson contribution, which is independent of £, is 
calculated using 



where S c = 250 mJy is the flux of the brightest source in our field of view. The uncertainty 

dN 
dS 



in the Poisson contribution involves the fourth moment J C S 4 jj£ dS of the source count 



and is given by 

2.32 



[AQ(0)]3 = I I 



/ q \ °- 84 

69.63 - 133.15 



(16) 



The analysis of large samples of nearby radio-galaxies has shown that the point sources 
are clustered. Cress et al. 1996 have measured the angular two point correlation function 
at 1.4 GHz (FIRST Radio Survey) , across an angular scale of .02° to 2°, equivalent to a 
£ range of 90 < £ < 9000. Throughout the entire angular scale the measured two point 
correlation function can be well fitted with a single power law of the form w(6) = (#/6 l o)~ /3 
, where P — 1.1 and 9 = 17 A'. This partly covers the range of angular scale (~ 10' to 20" 
or ~ 1000 < £ < 3 x 10 4 ) that we are interested in. We will assume that the clustering of 
the sources remain unchanged at our observing frequency. They have also reported that on 
small scales (< 0.2°) the double and multicomponent sources tend to have a larger clustering 
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amplitude than that of the whole sample. They also found that the sources with flux densities 
below 2mJy have a much shallower slope (~ 0.97) for the measured correlation function. It 
seems that the amplitude and slope of the measured two point correlation function changes 
with the angular scale and flux densities of the sources. For our present purpose we have 
used w(9) = (9/OqY 1 - 1 which have been measured up to I = 9000. We have assumed that 
the slope of the two point correlation function will remain unchanged beyond £ = 9000 . We 
then have 

where we oc £ /3_2 is the angular power spectrum which is the Fourier transform of w{6). 

The Galactic diffuse synchrotron radiation is believed to be produced by cosmic ray 
electrons propagating in the magnetic field of the Galaxy (Ginzburg & Syrovatskii 1969). 
The angular power spectrum is predicted to scale as £~ 7 with 7 ps 2.4 (Tegmark et al. 2000) 
to angular scales as small as 4' , and the spectral index has a value ~ 2.8. Here we have 
extrapolated the parameters from the 130 MHz model prediction of Santos et al. (2005). 
Recently Bernardi et al. (2009) have characterized the power spectrum of the total diffuse 
radiation at 150 MHz at the angular scales of our interest. The £ dependence that we adopt 
in our foreground model is consistent with that found by Bernardi et al. (2009). 

The total error in our model predictions is calculated by adding the variances from 
different contributions. 

Figure 7 shows the point source and synchrotron contributions along with the total mea- 
sured signal. At large angular scales (£ < 10 4 ) the foreground model prediction is dominated 
by the clustering of point sources, the point source Poisson contribution being the second 
largest component at these angular scales. This is reversed at smaller angular scales (£ > 10 4 ) 
where the point source Poisson contribution dominates and the clustering component is the 
second largest contribution. The Galactic synchrotron contribution, also shown in Figure 7, 
is much smaller at all the angular scales of our interest. The contributions from Galactic 
and extra-galactic free-free emission, whose parameters have been extrapolated from Santos 
et al. (2005), are also listed in Table 1. These are much smaller and hence are not shown 
in Figure 7. The expected HI signal {Ci ~ 10 -6 — 10~ 7 mK 2 ) is much smaller than all the 
foreground components mentioned here, and is not shown in the figure. 
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Figure 7. The black solid line shows the observed C'i with 1 — c error-bars, whereas the dash-dot-dot-dot curve shows the total 
foreground model prediction with its 1 — a error-bars. The dashed and the dot-dashed curves are respectively the Poisson and 
clustering contributions from the point sources, while the lower most dotted curve is the contribution from Galactic synchrotron 
radiation . 



We find that the measured Ct is within the 1 — a error-bars of the model prediction, 
for I < 2300. The measured Ct is around three times larger than the model predictions at 
smaller angular scales where the measured values do not lie within the 1 — a error-bars of 
the model predictions. The source of this discrepancy is, at present, unknown to us. The 
model predictions require the source properties to be extrapolated to faint flux levels and 
small angular scales where direct observations are not available. It is possible that the model 
predictions have been underestimated. For the present work we assume that the measured 
Ct is correct and that the model predictions have been underestimated at small angular 
scales. 

For the subsequent analysis in this paper we shall assume that the measured Ct(Au) is 
a combination of contributions from foregrounds, the HI signal and noise. Further, the HI 
signal being several orders of magnitude smaller than the foregrounds, we may interpret the 
measured Ct(Au) as an estimate of the foregrounds actually present in our field of view. 

6 FOREGROUND REMOVAL 

Removing the foregrounds which, as we have seen, are several orders of magnitude larger 
than the HI signal is possibly the biggest challenge for detecting the HI signal. There have 
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Table 1. Values of the parameters used for characterizing different foreground contributions at 610 MHz. Here S c is the flux 
of the brightest source in the field of view. 



Foregrounds 


yl(mK 2 ) 


a 


7 


Point source 
(Clustered part) 


/a „ \ 0.32 

20.03 x (%f ) 


2.07 


0.9 


Point source 
(Poisson part) 


8.38 x (Zjf-) 


2.07 





Galactic synchrotron 


0.122 


2.80 


2.4 


Galactic free-free 


1.14 x 10~ 4 


2.15 


3.0 


Extra Galactic free-free 


2.11 x 10~ 5 


2.1 


1.0 



been quite a few earlier works on this, nearly all either theoretical or simulation. All attempts 
in this direction are based on the assumption that the foregrounds are continuum radiation 
which vary slowly with frequency whereas the HI is a line emission which varies rapidly with 
frequency 

A possible line of approach is to represent the sky signal as an image cube where in 
addition to the two angular coordinates on the sky we have the frequency as the third di- 
mension. For each angular position, polynomial fitting is used to subtract out the component 
of the sky signal that varies slowly with frequency. The residual sky signal is expected to 
contain only the HI signal and noise (Jelic et al. 2008; Bowman, Morales & Hewitt 2009; 
Liu, Tegmark & Zaldarriaga 2009). Liu et al. (2009) show that this method of foreground 
removal has problems which could be particularly severe at large baselines if the uv sam- 
pling is sparse. They propose an alternate method where the frequency dependence of the 
visibility data is fitted with a polynomial and this is used to subtract out the slowly varying 
component. The residuals are expected to contain only noise and the HI signal. 

In this work we have attempted to subtract out the brightest point sources from the image 
using standard AIPS tasks. We have used the AIPS task 'UVSUB' to subtract the Clean 
Components (CC) of the brightest sources from the visibility data. Continuum images were 
used for this purpose. The resulting visibility data was used to make a new image. We find 
that this method fails to remove the point sources efficiently, several imaging artifacts remain 
in the vicinity of bright sources even after the sources have been removed. Similar findings 
were reported in Ali, Bharadwaj & Chengalur (2008) where the same technique was used to 
remove point sources from 150 MHz GMRT observations. Given the poor performance of this 
image based technique, we have not pursued it any further. The visibility based technique 
proposed by Liu et al. (2009) requires the data to be gridded in uv plane. The estimator 
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that we have used to determine Cg(Au) (Section 3) works with the individual visibilities. 
Using the gridded data would introduce a positive noise bias in Cg(Av) and hence we do 
not adopt this technique here. 

The foreground subtraction techniques discussed above all attempt to remove the fore- 
grounds before determining the angular power spectrum. Here we propose a different method 
where the foregrounds are subtracted after determining the angular power spectrum. The 
measured Cg(Av) (Figures 3 and 4) is a sum of the foregrounds, noise and the HI signal. 
The HI signal decays rapidly with increasing Av. This contribution is less than 10 % for 
Av > 0.5 MHz and it is negligibly small for Av > 1 MHz (Figure 6). We assume that 
Cg(Au) measured in the frequency interval 0.5 MHz < Av < 7.5 MHz contains only fore- 
ground and noise. Further, we assume that the foreground contribution to Cg(Av) has a slow 
Av dependence which can be well fitted by a low order polynomial. Note that, in addition 
to the intrinsic Au dependence of the foreground, the measured Cg(Av) has an additional 
Av dependence arising from the factor Q(Av) (eq. 4). The latter is a slow, monotonic vari- 
ation and we expect that both these effects can be adequately accounted for by a low order 
polynomial. We use the interval 0.5 MHz < Av < 7.5 MHz to estimate this polynomial, 
which is then used to subtract the foreground contribution from Cg(Au) across the entire 
range of our measurement (0 < Av < 7.5 MHz). The residual Ce(Av) is expected to be a 
sum of only the HI signal and noise. 

In order to illustrate our technique of foreground subtraction and to demonstrate its 
efficacy, we first apply it to simulated data where a known HI signal has been put in by 
hand. Given the uncertainty in our current understanding of the foreground properties and 
of the effects that have possibly been introduced during the observation and the subsequent 
analysis, we are guided by the measured Cg(Av) for our simulations. We find that the 
measured Cg(Av) (Figure 4) has a value around ~ 10 mK 2 , with ~ 5% variation with Av 
across the 7.5 MHz band. Further, the error has a typical value ^J[ACe(Av)] 2 ~ 0.01 mK 2 
(system noise only). We have simulated the measured MAPS using 

C e (Av) = £ a n (Av) n + 5 + a Cf\Av) (18) 

n 

where the polynomial J2n a n (Av) n represents the slowly varying Av dependence which 
causes Ce{Av) to vary by ~ 10 % across the 7.5 MHz band. Our Ce(Av) estimator (eq. 5) is 
even in Av, and hence we have only considered polynomials of even order. Our simulation was 
restricted to fourth order polynomials where the coefficients ao, 02, 04 are Gaussian random 
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Figure 8. The simulated Ci(Av) with 3<r error-bars (system noise) is shown for different values of I. The solid curve shows 
the best fit 4 th order polynomial determined using the interval 0.5 MHz < Au < 7.5 MHz. 



variables with mean 12, 0, mK 2 and rms. 1, 10 2 , 10 4 mK 2 respectively. The term 5 is a 
Gaussian random variable of rms. 0.01 mK 2 which incorporates the error and Q HI (Az/) is the 
HI signal (eq. 8) . The noise in our observation is considerably larger than the HI signal and it 
would not be possible to detect the signal even if the foregrounds were perfectly subtracted. 
The factor a in our simulations amplifies the HI signal so that it lies above the noise. The 
value of a has been chosen such that Cg(Ai/) = 5 x 0.01 mK 2 (5-sigma) at the value of 
Au where Ce(Au) is 70 % of the peak value Ce(0). The simulations have exactly the same 
frequency bandwidth and channel width as the measured Ci(Av). Though in this paper we 
have only considered fourth order polynomials for our simulations, the same procedure can 
easily be repeated considering even polynomials of any order. 
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Figure 8 shows the simulated Ce(Av) for the different values of £. Note that the poly- 
nomial coefficients a n are different for each realization of the simulation. We have fitted the 
simulated data with a fourth order polynomial using the interval 0.5 MHz < Av < 7.5 MHz. 
The best fit polynomial is also shown in Figure 8. The residuals, after the best fit poly- 
nomial is subtracted from the simulated Ci(Au), are shown in Figure 9. In the interval 
0.5 MHz < Av < 7.5 MHz, the residuals are within ±3a from zero which is consistent 
with noise. Figure 10 shows the residuals in the range Av < 1 MHz overlaid with the HI 
signal that had been added by hand. We find that our foreground subtraction technique 
successfully extracts the HI signal that had been added in the simulated data, despite its 
being buried in foregrounds which are ~ 200 times larger. We note that we have also tried 
a slightly different technique of foreground subtraction where we have used the entire Av 
range (< 7.5 MHz) to estimate the polynomial. We find that the latter technique does not 
correctly recover the HI signal that had been put in by hand. 

7 RESULT AND CONCLUSIONS 

We have measured the statistical properties of the background radiation across angular scales 
20" to 10' using the multi-frequency angular power spectrum Ce(Av). Frequency channels 
20 to 80 were used for the analysis. This corresponds to a total bandwidth of 7.5 MHz with a 
resolution of 125 kHz. The measured Ce(Av) has values around 12 mK 2 . Considering first the 
£ dependence of C t (Figure 3), starting from ~ 18 mK 2 at I ~ 1000, it drops to ~ 9mK 2 at 
I ~ 2000 and then rise to a nearly constant value of around 13 mK 2 . The uncertainty in C t is 
mainly due to the sample variance ie. the fact that we have observed a single ~ 1.5° x ~ 1.5° 
field of view which gives a limited number of independent estimates of Ce, the system noise 
makes a relatively smaller contribution. We next consider the Av dependence of Ce(Av) for 
different values of I (Figure 4). Assuming that the foreground contributions all have a smooth 
power law v dependence, the expected Av dependence may be estimated through a Taylor 



odd powers of Av/v cancel out because the estimator averages positive and negative Av 
values. The expected change in Ci(Av) is ~ 1.5 x 10~ 2 % for Av = 7.5 MHz. The measured 
Ci(Av) (Figure 4) has a smooth variation of the order of a few percent (1% to 4%) across the 
7.5 MHz bandwidth of our observation. In addition to the smooth Av dependence, we also 
notice a small oscillatory pattern in the measured Ce(Av). The expected HI contribution to 




where B is constant of order unity. The 
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Figure 9. The residual, with 3<r error bars, after subtracting the best fit 4th order polynomial from the simulated Ci(Au). 

Cg(Av) is ~ times smaller than the measured values, and we interpret the measured 
Ce(Av) as being nearly entirely foregrounds and noise. 

We next consider results for foreground removal using the technique discussed in Section 
6. For a fixed £, the frequency range 0.5 MHz < Av < 7.5 MHz was used to estimate a 
fourth order polynomial fit to Cg(Av). The Ci(-Av) = Ci{Av) symmetry of the Ci(Av) 
estimator was applied in the fitting procedure. This fit was used to subtract out the fore- 
ground contribution from the entire frequency range Av < 7.5 MHz. The performance of 
this foreground removal technique was assessed by visually inspecting the fit and the resid- 
uals across the entire band. We find that increasing the order of the polynomial does not 
result in any significant improvement, and hence we restrict our analysis to a fourth order 
polynomial for which the fits have been shown in Figure 4. The residuals in Ce(Av), we find, 
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Figure 10. The residual, with 3<r error bars, after subtracting the best fit polynomial from the simulated C'i(Au). The solid 
curve shows the HI signal that had been put in by hand in the simulations. 

typically have values within 0.1 mK 2 (Figure 11). In all cases the residuals are not consistent 
with Ci(Av) = (ie. noise only). The residuals, we find, have a nearly sinusoidal oscillatory 
pattern. These oscillations are most pronounced for the lowest £ value where it has an am- 
plitude of ~ 0.1 mK 2 and a period of Av ~ 3 MHz. The period and amplitude both decrease 
with increasing I. The oscillations are possibly not very well resolved at the larger £ values 
due to the 0.125 MHz channel width. The oscillations would possibly be more distinctly 
visible in observations with higher frequency resolution. The oscillatory residual pattern is 
quite distinct from the expected HI signal and also from random noise, and in principle it 
should be possible to distinguish between these by considering the Fourier transform 

C e (r m ) = Y J ^ Tm ^ Ct{Av n ) (19) 
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Figure 11. The residual, with 3c error bars (system noise only), after subtracting the best fit 4 th order polynomial from the 
measured C((Av). 

where n,m = -59, -58, 0, .., 58, 59, Au n = nxO.125 MHz and r m = m (119x0.125 MHz) -1 . 
We expect the oscillatory pattern to manifest itself as a localized feature in C^(r) and it 
should be possible to remove the oscillatory feature by applying a suitable filter to Q(r). 
We find that for the smallest £ values the amplitude of Ce(r) is peaked at a few r rn values 
located within I m \ < 10. Based on this we have chosen a filter 



F(r m ) = 
= 1.0 



-(H-m c ) 2 /2 



m \< m c 
m \> m r 



(20) 



such that F(r m )Ce(T m ) removes the Fourier components within | m |< m c from the residual 
Ce(r m ). Calculating Ci(Au) after applying the filter, we find that for the smallest £ value 
the oscillatory pattern is removed if we use m c = 7 or larger. The oscillatory pattern is 
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Figure 12. The residuals with 3<r error bars (system noise only) after applying the filter using m c = 7. 

somewhat reduced for the next two £ values while the three largest £ values are not much 
affected by the filter with m c = 7. It is possible to remove the oscillatory pattern from the 
second largest £ value by increasing the value of m c to m c = 14, but the oscillatory pattern 
still persists for the larger £ values. Increasing m c will also reduce the HI signal, and hence 
we do not consider m c = 14 in the subsequent discussion. The filter is also expected to affect 
the noise estimates, and the noise in the different Ci(Av) will be correlated as a consequence 
of the filter. For m c = 7, we are filtering out ~ 10% of the Q(r) values, and hence we do not 
expect this to be a very severe effect . Thus, for the purpose of this paper, it is reasonable 
to assume that the noise is unaffected by the filter. 

We find that for the smallest £ value (I = 1476) the residuals are consistent with zero at 
the 3a level. Based on this we conclude that we have successfully removed the foreground 
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Figure 13. The expected HI signal for xmb = 7.95 (circle and solid curve )and 3cr error bars (total errors) after applying the 
filter with m c = 7. The residuals from the measured C^(Af) are also shown (triangles). 



contribution from the measured Ce(Au) at this value of £. The residual oscillatory pattern 
persists at all the larger £ values where we are not successful in completely removing the 
foregrounds. The cause of this oscillatory residual, which at the moment is unknown to us, 
is an important issue which we plan to investigate in future. 

We next use the measured Ce(Au) at £ = 1476 to place an upper limit on the HI signal. 
The amplitude of the expected HI signal is determined by the factor (xmb) 2 (eqs. 8 and 
10) where xm and b are the HI neutral fraction and the HI bias parameter respectively. In 
the discussion till now we have used xmb = 2.45 x 10~ 2 to estimate the expected HI signal 
Cf l (Ai/). We now consider xmb as a free parameter whose value is unknown, and ask if it 
is possible to use our observation to place an upper limit on the value of xhi b. Considering 
xmb as an unknown parameter, the expected HI signal Cf l [xYub]{Au) can be expressed as 

* mb • ' (21) 



C?[x m b](Av) 



C?\Av) 



2.45 x 10- 2 

The HI signal would be detectable in our observation at the 3<j level if 



Cf l [x m b](Au) > ^{Cf l \x m b]{Au)Y/N E + {AC e (Au)} 



2 

sys 



(22) 



where Ne is the number of independent estimates of the signal, and the terms {C^^hi^] (Au)} 2 /Ne 
and {ACi(Av)} 2 sys are respectively the sample variance and system noise contributions to 
the total variance. 
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The fact that for I = 1476 the measured C^{Au) is consistent with noise, and the signal 
is not detected allows us to use eq. (22) to place an upper limit on Xmb- The filter F(t) 
that has been used to remove the oscillatory pattern in the residual also affects the signal. 
We have applied the same filter to Cf l [x m b}(Av) (Figure 13) and used this in eq. (22). The 
filtered signal is maximum at Au = and we use this data point to place a 3a upper limit 
on xmb. A value of xmb greater than 7.95 would have been detected in our observation, and 
is therefore ruled out at the 3a level. Our upper limit is around 330 times larger than the 
value that we have estimated based on results from quasar absorption spectra which imply 
Xm = 2.45 x 1CT 2 and the assumption that b = 1. The HI signal should, in principle, be 
detectable in observations that are a few hundred times more sensitive than the one that 
has been analyzed here. 
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